function grerr = hydro(filename)

  load(filename, 'params', 'tOut', 'uOut');
  params
  params.metricOrder = 1;


  % summary
  disp(sprintf('Grid Size   : (%d, %d, %d)', params.Nx, params.Ny, params.Nz));
  disp(sprintf('Grid Length : %f', params.Lx));
  disp(sprintf('Cell Length : %f', params.Dx));
  disp(sprintf('Target Time : %f', params.Tfinal));
  disp(sprintf('Time steps : %d',  params.Tsteps));
  disp(sprintf('Hydro Size  : %.2fMB', 3*params.Nx.*params.Ny.*5.*params.Tsteps/1024/1024));
  %disp(sprintf('Largest array  (X*X*Z*4*4*4) : %fMB', GRID_SIZE*GRID_SIZE*Z_SIZE*4*4*4/1024/1024));
  % switches





  % initialize derivative kernels
  h2dinit(params);




  

  Nx = params.Nx;
  Ny = params.Ny;
  Nt = Nx.*Ny
  size(tOut)
  size(uOut)

  params.Tsteps = 25;
  if true
    for t=1:params.Tsteps
      h2outplayer(reshape(tOut(t), 1,1), reshape(uOut(t,:), Nt*4, 1),'', params);
      u = squeeze(uOut(t,:));
      T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
      ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
      uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
      ut = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
      uxdx = dx(ux);
      uxdy = dy(ux);
      uydx = dx(uy);
      uydy = dy(uy);
      w = uxdy - uydx;
      max(max(abs(w)))
      %dlmwrite(sprintf('500-%d.txt', t), w, ' ');
      pause(1);
    end
  end

  return

  params.Tsteps = 25;
  %tOut = tOut(1000:1010,:)
  %uOut = uOut(1000:1010,:);

  if false
  % ****************************
  % compute time derivatives
  % ****************************

  tic;

  udtOut   = zeros(params.Tsteps, Nt.*5);
  udtdtOut = zeros(params.Tsteps, Nt.*5);

  % first derivative from eom
  for it=1:params.Tsteps
    t = tOut(it);
    u = uOut(it, :);
    
    udtOut(it,:) = h2eom(t, u, params);
  end

  % second der from finite differences, O(4)
  derker = [1.0/12; -2.0/3; 0; 2.0/3; -1.0/12];
  udtdtOut = conv2(udtOut, flipud(derker), 'same')./params.Dt;

  elapsed = toc;
  disp(sprintf('computed higher hydro time derivatives in %fsec', elapsed));

    global gr;


    tic
    perr = zeros(params.Tsteps, 4);
    for it=1:params.Tsteps
      u      = uOut(it,:);
      udt    = udtOut(it,:);
      udtdt  = udtdtOut(it,:);

      %perr(it, :) = h2constrel(u, udt, udtdt, params);
      % tic
      % h2metric(u, udt, udtdt, params);
      % metrictime = toc;
      % disp(sprintf('metric time step in %.3fsec', metrictime));
      % gr = -gr - 12;
      % mar1 = 2*params.Nz/32 + 1;
      % mar2 = 3*params.Nz/32 - 1;

      % integrate abs(gr) numerically 






      Nx = params.Nx;
      Ny = params.Ny;
      Nt = Nx*Ny;
      px = 3;
      py = 2;

      subplot(py, px, 3);
      plot(params.Cz, grz)
      title('-R-12');
      subplot(py, px, 4);
      plot(params.Cz, grz2)
      title('sqrt[g]*(-R-12)')

      subplot(py, px, 1);
      T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
      image((T-params.TempOffset)*params.colormapOffset+params.colormapOffset); 
      axis square;

      subplot(py, px, 2);
      image(reshape(gr(params.Nx/2,:,:), params.Ny, params.Nz)*params.colormapOffset+params.colormapOffset); 
      %axis square;







      refresh();
      pause(1);


    end
    toc
  end
  %figure;
  %plot(perr(:,1));


end


















